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AEROELASTIC  OPTIMIZATION  WITH 
MULTIPLE  CONSTRAINTS 

1.  INTRODUCTION 

This  report  summarizes  the  work  performed  under  the  sponsor¬ 
ship  of  the  Air  Force  Office  of  Scientific  Research,  with  Contract 
No.  F49620-78-C-0105.  The  overall  goal  of  this  research  has  been 
to  develop  improved  methods  of  designing  least-weight  aerospace 
structures  that  satisfy  both  strength  and  stiffness  requirements. 

During  the  first  year,  some  efficiency  comparisons  were 
undertaken  with  a  number  of  different  optimization  algorithms. 

All  of  these  algorithms  were  applied  to  problems  involving  a 
single  behavioral  constraint.  This  work  is  summarized  in  reference  1. 
From  this  work,  two  algorithms  were  chosen  for  work  with  multiple 
constraints.  The  first  algorithm  is  an  extension  of  work  performed 
by  Segenreich  and  co-authors  (refs.  2-4).  The  idea  here  is  to 
employ  slack  variables  in  order  to  transform  the  usual  inequality 
behavioral  constraints  into  equality  constraints.  In  this  manner, 
an  optimization  algorithm  based  on  multiple  equality  constraints 
can  be  applied.  Such  an  approach  has  been  applied  to  problems 
with  static  constraints  in  reference  4. 

The  second  algorithm  was  developed  by  Rizzi  (ref.  5) .  This 
algorithm  employs  a  Gauss-Seidel  iterative  procedure  to  aid  in 
the  elimination  of  inactive  constraints.  It  has  been  applied  by 
Rizzi  in  reference  5  to  problems  involving  both  strength  and 
stiffness  constraints,  and  it  was  felt  desirable  to  use  it  for 
comparison  with  the  first  algorithm  on  larger  problems.  Since 
both  of  the  algorithms  discussed  above  employ  a  similar  resizing 
formula,  the  comparison  proposed  is  primarily  a  test  of  two 
different  procedures  for  handling  multiple  constraints. 

The  sections  that  follow  will  describe  in  detail  the  two 
algorithms  as  coded  here,  the  numerical  applications  carried  out, 
and  some  recommendations  for  future  work. 


2.  STATEMENT  OF  PROBLEM  AND  THEORETICAL  BACKGROUND 


Consider  the  problem  of  minimizing  the  total  mass  M^  of 
a  structure.  Let  the  design  variables  be  t,f  i  =  l,2...,M, 
which  may  either  be  member-sizing  variables  in  a  finite-element 
model  of  the  structure  or  concentrated-mass  sizing  variables. 

With  the  usual  types  of  finite  elements,  the  mass  is  a  linear 
function  of  the  design  variables,  and  the  total  mass  can  be  written 
as 

Mt  -  Mo  ♦  Vi  (!) 


where  Mq  is  the  mass  not  available  for  optimization,  such  as 
fuel  mass.  There  are  behavioral  constraints  on  the  structure, 
representing  strength  or  stiffness  requirements  such  as  maximum 
allowable  stresses,  minimum  frequencies  of  free  vibration,  maximum 
displacement,  or  minimum  flutter  speeds.  These  are  characterized 
as  inequality  constraints 

cj(t±)  <0,  j  =  1,2, . . N  (2) 

Finally,  there  can  be  side  constraints  that  define  the 
maximum  and  minimum  values  of  the  design  variables: 


(t . )  .  <  t .  <  (t . )  , 

i  min  —  i  —  i  max 


i  =  1, 2, . . .  ,M 


(3) 


Necessary  conditions  for  a  local  minimum  are  given  by  the 
Kuhn-Tucker  conditions.  These  can  be  written  in  the  form 


N  3c. 
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fci  =  (ti}max 
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with  the  multipliers  Aj  given  by 


X3  =  ° 

if  Cj  <  0 

(7) 

Aj  >  ° 

if  Cj  =  0 

(8) 

There  are  essentially  two  parts  to  solving  an  optimization 
problem  as  posed  above.  The  first  involves  evaluating  the  constraints 
and  computing  derivatives  for  some  design?  the  second  involves  a 
redesign  algorithm  to  compute  the  next  design  change  in  an  iterative 
process  that  leads  to  a  local  optimum.  It  has  been  customary  in  the 
past  to  distinguish  between  direct  search  methods  (or  mathematical- 
programming  methods)  and  optimality-criterion  methods  for  redesign. 
However,  recent  work  (refs.  6  and  7)  has  shown  that  these  differences 
in  many  cases  are  more  philosophical  than  real.  For  present  purposes, 
it  is  sufficient  to  note  that  the  principal  problem  for  large  numbers 
of  constraints,  such  as  stress  constraints  on  individual  members,  is 
to  identify  in  an  efficient  manner  that  subset  of  the  constraints 
that  is  active  -  i.e.,  that  satisfies  equation  (8)  -  at  optimum. 
Subsequent  sections  will  describe  in  detail  two  strategies  that 
have  been  proposed  to  deal  with  this  problem. 

3.  AN  OPTIMALITY-CRITERION  METHOD  EMPLOYING 
SLACK  VARIABLES 


3.1  General  Description 

The  method  as  described  here  is  an  adaptation  of  that  described 
in  reference  4.  The  principal  idea  behind  this  method  is  to  use 
slack  variables  to  convert  inequality  constraints  to  equality 
constraints.  The  problem  of  identifying  active  and  inactive 
constraints  is  avoided,  and  the  redesign  algorithm  is  somewhat 
simplified,  since  it  need  be  set  up  only  for  equality  constraints. 

The  inequality  constraints  given  by  equation  (2)  are  written 
in  the  equivalent  form 


c !  (t . )  » 


cj (ti> 


+  z 
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These  constraints  are  appended  to  the  merit  function  (the  total 
mass)  in  the  form 


N 

M*  =  M  +  l  A.c! 
T  j=l  3  3 


=  + 


(10) 


Equations  (4)- (6)  remain  unchanged,  and  equations  (7)  and  (8) 
are  replaced  by 


AjZj  = 


j  - 


(ID 


c .  +  z?  =  0 
3  3 


/ 


3 


1,2, 


r  N 


(12) 


Equations  (11)  and  (12)  can  be  interpreted  to  mean  that  when  a 
constraint  c^  is  active  at  the  optimum,  then  the  corresponding 
slack  variable  z^  is  zero,  and  A^  can  have  an  arbitrary  value. 
Or,  if  the  constraint  is  inactive,  then  z..  is  not  zero  and  A^ 
must  be  zero  to  satisfy  equation  (11).  Equations  (4)- (6),  (11), 
and  (12)  constitute  2N  +  M  equations  for  2N+M  unknowns  -  the  M 
design  variables  t^,  the  N  multipliers  Aj  and  the  N  slack 
variables  z  ^ . 

The  redesign  algorithms  for  the  design  variables  and  the 
slack  variables  are  (ref.  4) 


f  v  V 

c:t;  , 

i  i 

(ti}min 

<  cvtv  <  (t.) 

—  ii—a.  max 

t?+1  = 
1 

< 

(ti}min' 

Citi  ^ 

lti>min 

(13) 

y,  Cti}maxf 

V  V 

c .  f*  > 

(ti>max 

zv+1  = 
3 

DV 

3 

(14) 
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where  v  is  the  iteration  number  and 


_v  v 
C .  =  a 

i 


+  (a  -1) 


(.15) 


dV  =  1  +  (aV-l) KX V  (16) 

The  redesign  factor  is  the  same  as  that  proposed  by  Kiusalaas 
(ref.  8).  From  the  optimality  criterion,  equation  (4),  it  can 
be  seen  that  cV  approaches  unity,  independently  of  the  value  of 
the  relaxation  parameter  a  ,  as  the  optimum  design  is  approached. 

From  equation  (11),  it  can  be  seen  that  the  factor  dV  approaches 
unity  in  similar  fashion  as  the  design  converges.  The  constant  K 
is  a  free  parameter  that  can  be  selected  by  the  designer;  in  effect, 
it  permits  controlling  the  step  size  on  the  slack  variables 
differently  from  the  step  size  on  the  design  variables. 

To  find  the  Lagrange  multipliers  X  ^ ,  the  requirement  that 
c !  =  0  for  all  j  is  imposed.  This  is  equivalent  to  the  procedure 
employed  in  an  optimization  problem  with  multiple  equality  constraints 
except  that  there  are  some  additional  terms  involving  the  slack 
variables.  Also,  for  the  sake  of  completeness,  terms  will  be 
included  to  account  for  design  variables  that  are  reaching  their 
upper  or  lower  bounds.  As  a  design  variable  reaches  one  of  these 
bounds,  it  is  left  at  the  appropriate  value  and  relegated  to  the 
passive  set.  After  a  converged  solution  has  been  obtained,  tests 
based  on  equations  (5)  and  (6)  are  applied  to  see  if  any  passive 
variable  is  to  be  reintroduced  to  the  active  set. 

Now  let  A^  denote  the  active  set  of  design  variables,  and 
assume  that  at  the  vth  iteration  equations  (13)  have  been  applied, 
and  those  design  variables  that  are  going  into  the  passive  set 
(not  those  relegated  to  the  passive  set  from  a  previous  iteration) 
have  been  identified.  Let  P^  and  P^  denote  such  design  variables 
for  the  lower  (or  minimum-value)  and  upper  (or  maximum-value) 
bounds,  respectively.  Then, 
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'  (cY-i) tY  ■ 

=  (ctV-l) 

r  ,  n 

l1  *  n  ih  x 

“I-  < 

(ti>»in  " 

‘i  • 

iEP«! 

(17) 

(t.) 

^  l  max 

*1  - 

iePv 

u 

To  first  order, 

(Ac ! ) V  =  AcV  +  2zVAzV  (18) 

3  3  3  3 

and 


4cj  -  •  iE(flv+P*+Pu) 


(19) 


After  use  of  equations  (17)  for  AtY,  equation  (19)  becomes 

*cj  -  +  l  elixl)  + 

A 

+  Jr  K’-*'1!]  (20) 


where 


(21) 


(22) 
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Equations  (14)  and  (16)  yield 


(ctV-l)KxVzV 
3  3 


(23) 


Equations  (20)  and  (23)  can  now  be  combined  to  give  an  expression 

for  (ACj)V.  Setting  this  quantity  to  zero  for  all  j,  letting 

(z?)v=-cY  ,  and  rearranging  the  resulting  equations  gives  a 
3  3  v 

linear  set  of  equations  for  the  Lagrange  multipliers  A ^ : 
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k=l 


‘3 
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'Sc .  \v' 

2v 

[(t . )  .  -t  .1 

L  i  mxn  ij 

/Sc .  \\j 

l  p  n 

+  .lv 
leP 
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[(ti)max-tIJ 
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(24) 


where 


Ki  - -  2Kc?kj  <25> 

and  6,  .  is  the  Kronecker  delta  function. 

*3 

Choosing  the  step-size  parameter  (aV-l)  completes  this 
algorithm.  There  are  four  criteria  that  are  applied?  these  follow, 
with  slight  modifications,  the  criteria  set  forth  in  reference  4. 
The  criteria  are  as  follows: 

(i)  Limit  on  step  size. 

The  change  in  each  design  variable  is  given  by 


,  v  , x  v^v 
(a  -l)viti, 


(26) 


where  v.  =  v./a^.  The  magnitude  of  the  vector  whose  components 

10  ii 

D 


are  AtY/tV  is 


(27) 


If  is 

K/tl  • 

where  NA 
Equation 
|av-l|: 


now  defined  as  the  maximum  permissible  magnitude  of 

1/2 

then  the  left-hand  side  of  equation  (27)  is  $t(NA)  7 
is  the  number  of  design  variables  in  the  active  set. 
(27)  can  then  be  rewritten  as  a  limiting  criterion  on 
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st(V1/2/ 


iv  nr] 

ieA  v  ' 

L  v  J 


1/2 


(28) 


(ii)  Feasibility  of  the  new  design. 

This  criterion  is  applied  in  order  to  ensure  that  a  design 
that  is  near  a  constraint  boundary  c^  =  0  will  not  be  driven  too 
far  from  it.  Equations  (12),  (18),  and  (23)  can  be  combined 
for  (Ac!)v  =  0  to  give  to  first  order 

Ac V  =  2 (aV-l) KxYcV  (29) 

With  6  defined  as  the  maximum  allowable  value  of  |AcV/cV|, 

C  3  J 

equation  (29)  becomes 


|aV-l|  1  6c/2K| xV |  (30) 

(iii)  Bounds  on  the  design  variables. 

This  criterion  is  used  to  ensure  that  no  design  variable  is 
driven  too  far  outside  its  allowable  range.  Tolerance  bands  on 
the  upper  and  lower  bounds  are  established  with  parameters  t^(>l) 
and  t^(<l),  respectively.  With  the  proviso  that  (av-l)  be  always 
negative,  this  criterion  takes  two  forms.  From  equation  (26) , 
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it  can  be  seen  that  the  upper  or  lower  bound  applies  if  vY 
is  respectively  negative  or  positive.  Hence,  the  criterion  is 
one  of  the  following: 


<ti>min]/viti 

(ti,max]/vlti 


vY  >  0 
1 

vv  <  0 

l 


(31) 

(32) 


( iv )  Convergence 

As  the  design  nears  an  optimum,  it  is  possible  that  the 
criteria  outlined  above  will  all  provide  relatively  large  upper 
bounds  on  |aV-l|.  To  ensure  that  the  algorithm  remains  stable, 
this  criterion  merely  establishes  an  upper  limit  on  |av-l|: 


av-l|  <  d 


(33) 


At  any  step  in  the  iterative  process,  the  value  of  (a  -1) 
to  be  used  is  given  by 


(34) 


where  the  right  side  is  given  by  the  minimum  of  the  upper  bounds 
computed  from  equations  (28),  (30),  (31),  (32),  or  (33). 

There  are  additional  details  that  need  to  be  discussed 
concerning  the  assignment  of  a  design  variable  to  the  passive 
set.  The  parameters  t^  and  ty  are  used  to  define  tolerance  bands 
on  both  upper  and  lower  bounds.  In  other  words,  a  design  variable 
is  consigned  to  the  passive  set  and  given  the  appropriate  upper- 
or  lower-bound  value  whenever 


(36) 
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or 


v+1 


t.  VCt.)  >  t0 
x  '  1  max  —  i 


(37) 


It  is  intended  here  to  accelerate  as  much  as  possible  the  selec¬ 
tion  of  active  and  passive  variables.  When  t  and  t0  are  not 
equal  to  unity,  this  means  that  a  design  variable  may  be  selected 
as  passive,  even  though  (ct  -1)  has  not  been  selected  from  equa¬ 
tions  (31  or  (32).  As  a  consequence  of  this,  either  of  the 
sets  P^  or  Pv  will  no  longer  be  null,  and  it  will  be  necessary 
to  return  to  equation  (24)  and  repeat  the  calculation  of  the 
Lagrange  multipliers.  In  principle,  this  process  must  be  repeated 
until  the  selection  of  the  passive  variables  converges;  in 
practice,  only  one  repetition  is  usually  necessary.  Note  also 
that  if  tf  and  t  are  set  equal  to  unity,  then  the  choice  of 
(a  -1)  will  be  governed  by  one  of  the  equations  (31)  or  (32)  and 
no  repetition  will  be  needed,  since  there  is  no  inconsistency  in 
equations  (17),  which  give  the  increments  AtV. 

Since  all  of  the  a.'s  are  positive,  convergence  to  a  local 
minimum  is  measured  by  the  v^,  and  the  convergence  criterion  is 

|vY|  1  ev  '  ieAv  (38) 


where  ev  is  a  user-supplied  quantity.  To  ensure  that  no  passive 

design  variables  need  to  be  reintroduced  into  the  active  set,  the 

applicable  Kuhn-Tucker  conditions,  equations  (5)  and  (6),  are 

verified  (with  v?  rather  than  vY): 

l  i 

v^  >  0  ,  ieP^  (39) 


0  , 


(40) 
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Here  and  P^  define  the  complete  set  of  passive  design  variables. 
Since  the  redesign  formula  gives  equation  (26)  for  At?  ,  there 
is  an  obvious  physical  explanation  for  the  above  conditions.  If 
a  design  variable  t.  is  at  its  minimum  value  and  its  associated 
convergence  function  v^  is  negative,  a  redesign  will  produce  a 
positive  increment  At^,  thereby  implying  that  this  design  variable 
should  be  reintroduced  into  the  active  set.  Similar  reasoning 
holds  for  a  design  variable  at  its  upper  bound  when  its  conver¬ 
gence  function  is  positive. 


3.2  Summary 

The  algorithm  is  summarized  in  the  following  steps; 

(1)  For  the  current  design,  evaluate  the  constraints  c . ,  j=l,2,...,N. 

(2)  Calculate  the  derivatives  3Cj/3ti  ,  i  =  l,2,...,M;  j=l,2,...,N. 

(3)  Calculate  3^  from  equation  (21)  and  F.  •  from  equations  (22) 
and  (25). 

(4)  Solve  for  the  from  equation  (24),  with  P^  and  Pu  null 

if  this  is  the  first  time  the  multipliers  are  calculated  for 
this  iteration.  If  not,  include  terms  for  the  appropriate 
design  variables  identified  in  step  (8)  below. 

(5)  Compute  the  convergence  functions  v^  =  v\/a^  '  vi  9-i-ven  by 
equation  (4) . 

(6)  Test  for  convergence,  using  equation  (38).  If  this  test  is 
satisfied,  test  the  passive  variables  with  equations  (39)  and 
(40).  If  any  of  these  tests  fails,  redefine  the  active  and 
passive  sets  appropriately  and  return  to  step  (3) .  If  the 
tests  are  satisfied,  exit  with  the  final  design  information. 

(7)  Compute  (av-l)  from  equations  (28)  and  (30)- (34). 

(8)  Compute  the  new  design  variables  according  to  the  first  of 
equations  (13)  and  equation  (15) .  Test  for  new  passive 
variables  using  equations  (36)  or  (37).  If  necessary,  redefine 
the  active  and  passive  sets  and  return  to  step  (3) .  Other¬ 
wise,  compute  the  new  mass  and  return  to  step  (1)  to  begin  a 
new  iteration. 
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In  brief,  this  algorithm  seeks  to  avoid  identifying  active 
constraints  by  using  slack  variables  to  permit  all  of  the  constraints 
to  be  treated  as  equality  constraints.  Implied  in  this  approach 
is  the  assumption  that  all  intermediate  designs  and  the  starting 
design  must  be  feasible  -  that  is,  no  constraint  violations  are 
tolerated.  Also,  all  constraints  are  to  be  evaluated  at  every 
design,  as  the  algorithm  is  currently  formulated.  For  truly  large 
systems,  where  a  great  number  of  constraints  may  be  applied,  it 
would  be  necessary  to  introduce  a  "throwaway"  concept  or  some 
similar  strategem,  whereby  constraints  that  are  nowhere  near  being 
active  are  at  least  temporarily  ignored. 

4.  AN  OPTIMALITY-CRITERION  METHOD  EMPLOYING 
GAUSS-SEIDEL  ITERATION 

4.1  General  Description 

This  method  is  based  on  the  work  of  reference  5.  The  stra¬ 
tegy  here  is  to  let  the  active  constraints  be  selected  by  the 
signs  of  the  Lagrange  multipliers.  The  Gauss-Seidel  iterative 
procedure  is  used  to  permit  the  constraints  associated  with 
negative  multipliers  to  be  eliminated  efficiently.  At  the  optimum, 
the  only  constraints  involved  are  the  active  ones,  for  which  the 
Lagrange  multipliers  are  non-negative. 

The  problem  statement  is  given  by  equations  (l)-(8),  with 
the  redesign  algorithm  given  by  equations  (13)  and  (15) .  In  order 
to  reduce  the  number  of  constraints  included  in  the  calculation 
of  the  A-,  a  positive  tolerance  parameter  6^  is  specified,  and 
those  constraints  are  selected  for  which 


Let  A^  define  the  set  of  constraints  for  which  this  test  is  satis- 
c 

fied.  These  constraints  are  either  potentially  active  or  violated. 


Computation  of  the  first-order  approximation  for  AcV  gives 
equation  (20).  If  it  is  now  assumed  that  all  the  constraints  in 
are  to  be  active  constraints,  then  it  is  required  that 
AcV  =  -Cj  ,  and  equation  (20)  can  be  manipulated  to  give 


~  '  °i  +  UP?  (^) 


&  (%)v 


/  (av-l) ,  jtA^ 


The  coefficients  fiV  and  8^.  are  given  by  equations  (21)  and  (22) 
respectively,  with  j,keAc. 

The  next  step  is  to  determine  that  subset  of  AV  for  which 

c 

solutions  of  equations  (42)  will  produce  non-negative  values  of 
X^.  One  approach  is  to  solve  equations  (42)  for  the  set  of  A^, 
eliminate  the  terms  associated  with  any  negative  multipliers, 
and  repeat  this  process  until  only  non-negative  values  of  the 
multipliers  are  computed.  However,  as  noted  in  reference  8,  there 
is  no  guarantee  that  this  procedure  will  converge,  because  of  the 
coupling  among  the  equations  provided  by  the  coefficients  8^j  , 
k  /  j.  The  solution  proposed  in  reference  (5)  is  to  apply  the 
Gauss-Seidel  iterative  procedure  to  equations  (42).  Let  these 
equations  be  rewritten  in  matrix  form  as 


[B] { A )  =  {b> 


(43) 


and  let  n  be  their  order.  Then,  let  equation  (43)  be  recast  in 
the  following  form: 


(X)  =  [B]{X)  +  (b) 


(44) 


where 


B^j  —  —  ^  l#2,...,n  (45) 

=  b^/B^  ,  i  =  lf2fa«»fn  (46) 

Thus,  [B]  is  a  matrix  with  a  zero  principal  diagonal,  and  the 
Gauss-Seidel  iterative  algorithm  is 

V+  V  B. .xV-1  +  b.  ,  i  =  1, 2, . .  .  ,n  (47) 

D  j —  i+1  «  3  i 

Here  y  is  the  Gauss-Seidel  iteration  number,  and  the  terms  in  the 
first  summation  for  i  =  1  and  the  second  summation  of  i  =  n  are 
ignored.  As  this  algorithm  is  applied,  each  multiplier  X^  is 
estimated  from  current  values  for  i  <  k.  This  means  that  when¬ 
ever  a  negative  value  for  a  multiplier  is  determined,  it  can  be 
immediately  deleted,  with  minimal  impact  on  the  values  computed 
for  the  remaining  multipliers  in  that  particular  iteration. 
Although  it  cannot  be  rigorously  proven,  this  process  is  expected 
to  lead  naturally  to  the  proper  subset  of  A  (ref.  5).  Let  this 

subset  be  AJ.. 

c 

Redesign  can  now  be  carried  out  with  equations  (13)  and  (15) 

with  jeAv  in  equation  (15) .  The  process  of  identifying  passive 
c 

design  variables  is  similar  to  that  described  in  Section  3  above. 
If  any  design  variables  are  predicted  to  reach  their  upper  or 
lower  bounds,  then  the  calculations  of  the  multipliers  must  be 
repeated,  with  the  appropriate  terms  added  to  equation  (42) . 
Unlike  the  original  algorithm  proposed  in  reference  5,  design 
variables  remain  passive  until  a  converged  design  is  found.  At 
that  point  the  tests  given  by  equations  (39)  and  (40)  are  applied 
and  the  appropriate  variables  reintroduced  into  the  active  set  if 
any  of  the  tests  fails. 


To  accelerate  convergence,  a  can  be  chosen  as  a  positive  number 

(OT 

greater  than  unity,  with  a  given  as  a  negative  number.  The 
multiplier  6X  is  typically  less  than  unity,  since  it  is  desirable 
to  have  equation  (41)  provide  an  increasingly  restrictive  test  as 
the  optimization  proceeds. 

Convergence  is  measured  by  two  criteria: 


|Cj|  <  ec  ,  jeA^  (50) 


< 


(51) 


4 . 2  Summary 

This  algorithm  is  summarized  as  follows: 

Cl)  For  the  current  design,  evaluate  the  constraints  c j , 

3  =  1, 2, .  . .  ,N. 

(2)  Calculate  the  derivatives  3c. /3t.,  i  =  l,2,...,M,  J-1,2,...,N. 

J  ^  \) 

(3)  Use  equation  (41)  to  test  for  the  set  A  of  potentially  active 
or  violated  constraints. 

(4)  Compute  and  p^  .  from  equations  (21)  and  (22),  respectively, 
with  j,keA^. 

(5)  Use  the  Gauss-Seidel  iterative  formula,  equation  (47)  to  deter¬ 
mine  the  constraint  set  A^. 

(6)  Compute  the  redesign  factors  cV  from  equation  (15)  and  identify 
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any  new  passive  variables,  redefine  the  active  and  passive 
sets  accordingly  and  return  to  step  (4),  with  the  appropriate 
terms  added  in  equation  (42).  Continue  this  loop  through 
steps  (4)- (6)  until  there  is  no  change  in  the  active-passive 
identities. 

(7)  Compute  the  convergence  functions  vV  for  the  constraint  set 

Av  and  test  for  convergence  with  equations  (50)  and  (51) . 

c 

If  these  tests  are  satisfied,  test  the  passive  variables  with 
equations  (39)  and  (40).  If  any  of  these  tests  fails,  redefine 
the  active  and  passive  sets  accordingly  and  return  to  step  (4). 
If  not,  exit  with  the  final  design  information. 

(8)  If  any  of  the  convergence  tests  fails,  compute  the  new  design 

variables  (equations  (13))  and  the  new  values  of  the  step-size 

parameter  a  and  the  constraint-tolerance  parameter  6” 

c 

(equations  (48)  and  (49))  and  return  to  step  (1)  to  begin  a 
new  iteration. 

It  is  worth  noting  that  this  algorithm  also  requires  that  all 
constraints  be  evaluated  at  each  iteration.  The  selection  of  the 
active  constraints  -  that  is,  those  that  are  actually  included  in 
the  redesign  -  is  a  two-stip  process.  First,  the  constraints  that 
are  violated  or  are  within  a  given  tolerance  band  around  zero  are 
chosen,  and  from  these  the  active  set  is  chosen  by  eliminating 
those  for  which  the  multipliers  are  negative.  It  is  therefore  to 
be  expected  that  this  strategy  would  be  most  valuable  when  there 
is  a  large  number  of  constraints,  and  the  active  set  is  rapidly 
changing  from  iteration  to  iteration.  When  the  constraints  are 
so  numerous  that  it  is  not  feasible  to  evaluate  them  all  at  every 
iteration,  a  "throwaway"  concept  would  be  useful  and  could  easily 
be  incorporated. 

5.  CONSTRAINT  EVALUATIONS  AND  DERIVATIVE  CALCULATIONS 

Three  types  of  constraints  are  considered  -  displacement  in 
response  to  a  fixed  static  load,  natural  frequencies,  and  flutter 
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speed.  The  evaluations  of  the  constraints  and  their  derivatives 
follow  generally  the  procedures  outlined  in  reference  1.  The 
basic  formulas  will  be  repeated  below  for  the  sake  of  completeness 
and  any  differences  between  current  usage  and  that  in  reference  1 
will  be  described. 


5.1  Displacement  Constraint 
The  constraint  is  given  by 


c 


u 


<ur>des 


1 


(52) 


and  the  derivative  is  calculated  by  the  dummy-load  method: 

,  3u 


3c 


lur>des 


r 
3 1 . 


tur>des 


[uj  [K.  ]  (u(r)} 


(53) 


Here  (u)  is  the  set  of  nodal  displacements  associated  with  a  given 
loading  condition,  whose  rth  component  u  is  constrained  to  be  less 
than  or  equal  to  (ur)des  »  and  {uv  is  the  set  of  displacements 
resulting  from  the  application  of  a  dummy  load  at  the  rth  node  in 
the  direction  of  ur.  The  matrix  [K^]  is  the  stiffness  matrix 
associated  with  design  variable  t.;  both  mass  and  stiffness 
matrices  vary  linearly  with  the  design  variables.  For  multiple 
constraints,  different  loading  conditions,  and  different  displace¬ 
ments  for  each  loading  condition,  can  be  considered. 


5.2  Frequency  Constraints 

In  contrast  to  the  expression  in  reference  1,  the  constraint 
function  for  a  natural-frequency  constraint  is  written  as 


W 


'“r'des 
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(54) 


where  w  is  the  rth  free-vibration  frequency  of  the  structure,  and 
it  is  constrained  to  be  greater  than  or  equal  to  (a  . 

The  derivative  expression  is 


^°r^  des 


9  (o 

_ r 

5T7 

1 


(55) 


where 


lq(r)  J  [GKi]{qtr)>  -  Lq  Cr  ]  J  [GMi]  {q  (r)  } 
2<o  lq  ^r)  ]  [GM]  (q  {r)  } 


(54) 


It  is  assumed  here  that  the  natural  modes  of  the  initial  design 
are  used  to  define  generalized  coordinates  for  all  the  designs,  so 
that  the  generalized  derivative  matrices  (GK^l  and  [GM^l  can  be 
written  as  an  inactive  part  plus  the  sum  of  the  active  parts. 

The  column  {q  '}  is  the  rth  mode  shape  of  the  current  design 

(see  reference  1).  Multiple  frequency  constraints  can  be  considered. 

5.3  Flutter-Speed  Constraints 

Here,  the  constraint  function  c  is  identified  directly  with 
the  damping  parameter  g.  In  essence,  this  formulation  requires 
that  the  imaginary  part  of  a  critical  root  of  the  flutter  equations 
be  non-positive  at  a  given  free-stream  speed  and  altitude.  The 
mode  shapes  of  the  initial  design  are  used  to  define  generalized 
coordinates,  so  the  generalized  derivative  matrices  (GK^]  and 
(GM^l  are  constant,  and  the  generalized  aerodynamic  matrix  [GA] 
varies  only  with  reduced  frequency  k.  The  derivative  is  evaluated 
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just  as  in  reference  1: 


8c  8g 

8 1 .  8 1 . 
l  l 


r>i  2_i 

R2  “  u  R^ 
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3  +  2I3  + 


3 

-  (2R3  -  2gl3  +  r4)(I2  “  “2ll  +  9R 2 )  /D  (57) 


where 


=  (2gR,  +  21,  + 


k^L.  1  \ 
v-  4/ 


I3  +  ^2R3  -  2gl3  +  *2-  R4JR3  (58) 


Rj  -  ReflpJlGMjHq})  (59) 

R3  =  ReflpJ  (GKi]{q})  (60) 

R3  =  Re  (ipj tGK]{q})  (61) 

r4  -  Re(lpj[^](i})  (62) 

and  1^,  I^f  I3,  and  I4  are  corresponding  imaginary  parts.  Also, 
is  the  free-stream  speed,  b  is  a  reference  length,  a  is  the 
frequency,  {q}  is  the  eigenvector  for  the  critical  flutter  mode, 
and  {p}  is  the  adjoint  eigenvector. 

When  the  constraint  is  evaluated,  it  is  necessary  to  ensure 

that  the  reduced  frequency,  which  is  fixed  in  order  to  evaluate 

[GA] ,  is  compatible  with  the  frequency  computed  from  the  eigenvalue 
—  2  1/2 

ft  =  (l+ig)/w  (here  i =  (-1)  ) .  While  it  is  possible  to  compute 

and  use  it  to  estimate  the  new  value  of  the  frequency  once  the 
design  change  has  been  compuated,  it  is  simpler  to  iterate.  The 
previous  value  of  reduced  frequency  is  used  to  define  the  aero- 
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dynamic  matrix  [GA] ,  and  a  value  of  w  computed  from  the  eigenvalue 
TT.  If  this  value  is  not  equal  to  within  some  specified  tolerance 
to  the  value  used  to  compute  the  reduced  frequency,  it  is  used  as 
the  new  reduced  frequency  and  the  matrix  equation  for  flutter  is 
solved  again.  Typically  only  a  few  iterations  are  needed;  more 
details  can  be  found  in  reference  1. 

Multiple  flutter-speed  constraints  may  be  considered.  These 
may  be  critical  speeds  at  different  altitudes  or  Mach  numbers,  or 
they  may  involve  additional  potentially  critical  branches  for  a 
single  case. 


5.4  Accuracy  of  Derivative  Calculations 

In  order  to  test  the  accuracy  of  the  derivative  calculations, 
the  analytical  derivatives,  computed  as  described  above  for 
frequency  and  flutter  constraints,  were  compared  with  derivatives 
computed  by  finite  differencing.  The  structure  used  was  the 
sandwich  wing  of  reference  9.  This  wing  had  a  primary  flutter 
constraint  corresponding  to  M^^O.8  at  an  altitude  of  2,438  m, 
a  secondary  flutter  constraint  corresponding  to  the  crossover  of 
another  root  in  the  V-g  plane  at  a  slightly  higher  airspeed,  and 
a  frequency  constraint  that  required  the  fundamental  natural 
frequency  to  be  not  less  than  that  of  the  initial  design.  Addi¬ 
tional  details  concerning  this  wing  and  the  constraints  are  given 
in  Section  6. 

Finite-difference  derivatives  with  respect  to  all  seven  design 
variables  are  compared  to  analytical  derivatives  in  Table  1.  These 
derivatives  were  computed  at  an  intermediate  design  obtained  from 
an  optimization  run,  so  the  modes  used  to  define  generalized  coor¬ 
dinates  are  primitive  modes.  As  can  be  seen,  the  derivatives 
calculated  by  the  two  methods  are  generally  in  very  good  agreement. 
Finite-difference  derivatives  were  computed  for  a  number  of  different 
perturbation  magnitudes  and  compared  in  order  to  minimize  numerical 
errors;  the  derivatives  given  in  the  tables  were  determined  with 
5%  perturbations. 
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During  this  portion  of  the  study,  it  was  discovered  that  the 
convergence  test  used  in  the  iterating  to  determine  the  new  frequency 
when  the  flutter  constraint  was  being  evaluated  had  to  be  tightened 
considerably  in  order  to  ensure  sufficiently  accurate  determination 
of  the  damping  parameter  g.  It  was  ultimately  necessary  to  require 
that  the  selected  and  computed  values  of  reduced  frequency  agree 
to  within  0.0001. 

As  noted  in  reference  10  and  elsewhere,  the  use  of  fixed  modes 
does  lead  to  inaccuracies  both  in  constraint  evaluation  and  in 
derivative  calculation.  On  the  other  hand,  updating  the  modes 
after  every  iteration,  or  even  after  a  certain  number  of  iterations, 
is  often  not  feasible  because  of  the  additional  computer  time 
required  or  the  abrupt  changes  introduced  in  constraints  and 
derivatives  when  the  modes  are  changed.  In  principle,  using  a 
sufficiently  large  number  of  modes,  even  if  they  are  used  as 
fixed  modes,  should  resolve  this  difficulty.  However,  there  are 
practical  limits  to  the  number  of  modes  that  can  be  used.  For 
example,  mode  shapes  are  often  fit  with  polynomials  or  splines 
before  they  are  used  in  calculating  generalized  aerodynamic  forces, 
and  the  accuracy  of  such  fits  for  high-order  modes  with  complicated 
node-line  patterns  is  questionable. 

Another  approach  has  been  proposed  in  reference  5.  It  was 
hypothesized  that  the  principal  source  of  modeling  errors  with 
fixed  modes  is  the  failure  of  the  inplane  portions  of  the  modes 
to  produce  negligible  inplane  forces  as  the  design  is  changed. 

To  avoid  this  problem,  it  was  proposed  to  preserve  unchanged  the 
transverse  displacements  but  recompute  the  inplane  displacements 
from  the  requirement  that  the  inplane  forces  be  zero.  This  would 
involve  additional  computations  in  constructing  the  generalized 
mass  and  stiffness  matrices  but  would  not  affect  the  generalized 
aerodynamic  forces.  In  reference  5,  the  inplane  displacements 
were  recomputed,  but  the  additional  terms  that  arise  in  the 
derivative  calculations  to  account  for  changes  in  certain  portions 
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of  the  modes  were  ignored.  The  accuracy  of  constraint  evaluation 
was  enhanced  by  this  strategem,  and  it  is  expected  that  accounting 
for  the  mode-shape  changes  in  the  derivative  calculations,  which 
involves  only  the  derivatives  of  the  generalized  mass  and  stiffness 
matrices,  should  improve  the  accuracy  of  the  derivative  calculations 
with  little  added  computational  effort. 

6 .  NUMERICAL  EXAMPLES 

6.1  Test  Cases 

It  is  instructive  to  consider  first  some  test  cases  involving 
two  design  variables,  so  that  the  progress  of  the  optimization  can 
be  plotted  and  visualized  in  design  space.  For  all  of  the  test 
cases,  the  merit  function  was 

M  =  2tx  +  3t2  (63) 

For  the  first  test  case,  the  "behavioral"  constraints  were 

C1  =  t2  ~  fcl  -  0  (64) 

c2  -1  -  (t*  +  t2  )/2  <  0  (65) 

and  minimum-gage  constraints  of  (t. )  .  =0.2,  (t0)  .  =0.5  were 

also  imposed.  Figure  1  presents  the  path  in  design  space  taken 
by  the  slack-variable  program  (program  SLACK)  from  an  initial 
design  at  (2. 0,1. 5).  As  can  be  seen  the  first  few  steps  follow 
reasonably  closely  a  steepest-descent  path,  nearly  normal  to  the 
M  =  const,  lines,  until  the  minimum-gage  constraint  on  t2  is 
encountered.  Thereafter  the  designs  proceed  along  this  constraint 
boundary  until  the  intersection  with  the  c2  =  0  boundary  is  reached, 
at  (1.3228,  0.5000).  The  value  of  the  merit  function  at  this  pointf 


which  also  is  a  global  minimum,  is  (M)  .  =  4 . 1456 .  After  42  itera- 

tions/  the  computed  values  for  the  optimum  were  CM)Jn^n  =  4,14  60  at 
(1.3230,  0.5000).  Values  of  other  parameters  used  in  SLACK  are  given 
in  the  caption.  No  attempt  was  made  at  this  point  to  vary  these 
parameters  so  as  to  reduce  the  number  of  iterations  required. 

The  next  set  of  slack  test  cases  was  run  with  the  sign  on  c^ 
changed,  so  that  the  feasible  domain  was  to  the  left  of  the  c^  -  0 
boundary.  With  the  same  minimum-gage  constraints,  there  are  multiple 
extrema  here  -  one  at  (0.2000,  1.4000),  with  = 4 . 600,  and  one 

at  (1.0000,  1.0000),  with  (M)  = 5 .000 .  With  the  initial  design  at 

(0.9000,  1.2000),  the  latter  extremum  is  found  after  44  iterations 
(figure  2);  the  computed  values  are  (M)  = 5. 0019  at  (0.9981,  1.0020). 

If  the  initial  design  is  moved  to  (1.2000,  2.000),  the  other  extremum 
is  found  after  64  iterations  (figure  3).  The  computed  values  here 
are  (M)  .  =4.6006  at  (0.2000,  1.4002).  In  figures  2  and  3,  a  line 

is  drawn  through  the  origin  in  the  steepest-descent  direction.  The 
intersection  of  this  line  with  the  c2  =  0  boundary  (point  A)  defines 
the  maximum  value  of  M  on  the  boundary.  It  is  therefore  readily 
apparent  that  any  optimization  step  reaching  the  c2  =  0  boundary  to 
the  right  of  point  A  will  result  in  an  optimum  solution  to  the  right, 
and  in  similar  fashion  the  optimum  to  the  left  will  be  found  for  any 
step  reaching  the  c2  =  0  boundary  to  the  left. 

Another  test  case  with  the  same  constraints  and  c2,  but  with 

(t..)  .  =  0.500,  (t0)  .  =1.500,  was  run  to  test  the  operation  of 

1  min  i  mm 

SLACK  when  the  optimum  design  is  dictated  by  the  minimum-gage 
constraints.  This  optimum  was  found  after  17  iterations  (figure  4) . 

The  final  SLACK  test  case  involves  two  concave  constraints  - 
the  same  constraint  function  c1  used  above,  and 

c2  =  (tx-2)2  +  (t2-2)2  -  2.56  <  0  (66) 

Minimum-gage  constraints  of  (0.4000,  0.4000)  were  imposed,  but 
these  were  inactive,  and  the  global  optimum  is  at  (0.9102,  0.8285), 


with  (M)  .  *  4.3060  Cfigure  5}.  With  an  initial  design  of 

mm 

(1.000,  2.000),  the  optimum  was  found  after  32  iterations.  The 

computed  values  were  (M)  =  4.3084  at  (0.9072,  0.8313), 

min 

It  is  interesting  to  note  that  the  behavior  of  the  slack- 
variable  algorithm,  as  given  by  the  paths  in  design  space  plotted 
for  these  test  cases,  closely  follows  a  gradient-projection 
strategy:  steepest  descent  until  a  constraint  surface  is  encountered, 
and  then  steps  along  the  projection  of  the  gradient  of  the  merit 
function  on  the  tangent  to  the  constraint  surface  (see,  for  example, 
reference  11).  This  is  not  so  surprising,  in  view  of  the  close 
relationship  between  a  Kiusalaas-type  algorithm  and  a  projection 
algorithm  (reference  12). 

A  number  of  cases  with  different  parameters  in  SLACK  were 

run  with  the  constraints  and  initial  design  given  in  figure  5. 

The  results  are  summarized  in  Table  2.  As  would  be  expected, 

imposing  a  more  severe  convergence  criterion  results  in  an 

increased  number  of  iterations  required  -  in  this  case,  the  number 

of  iterations  is  almost  doubled  for  a  factor  of  10  reduction  in 

ev.  Comparing  numbers  of  iterations  with  different  values  of  K 

and  6  suggests  that  larger  values  of  6  and  smaller  values  of 
c  c . 

K  will  enhance  convergence,  at  least  in  cases  where  many  iterations 
take  place  along  a  nonlinear  constraint  boundary.  This  is  in  fact 
expected  to  occur  for  most  problems  of  practical  significance. 

One  test  case  was  run  with  the  program  based  on  Gauss-Seidel 
iteration  -  program  GAUSS.  This  case  involves  the  same  merit 
function  as  previously  used  and  three  constraints: 

C1  =  fcl  “  t2  I  0  (67) 


c2  =  -1  +  (t^  +  t2)/2  <  0  (68) 

c3  =  2  (l”t^  )  -  t2  <  0  (69) 
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These  constraint  boundaries  and  the  results  from  a  run  with  program 
GAUSS  are  presented  in  Figure  6.  The  feasible  domain  is  a  closed 
area  roughly  triangular  in  shape,  with  a  global  optimum  at  the 
intersection  of  the  boundaries  of  the  two  active  constraints,  c^ 

and  c,.  The  optimum  is  at  (0.8165,  0.6667),  with  (Ml  .  = 3. 6330. 

J  nun 

With  an  initial  design  at  (0.8200,  0.9300),  the  optimal  solution 

was  found  in  just  four  interations  to  be  (M)  =3.6718  at 

min 

(0.8115,  0.6830).  As  can  be  seen  in  the  figure,  there  were  some 
oscillations  about  the  optimum.  These  were  undoubtedly  the  result 
of  the  values  used  for  the  various  parameters,  which  are  given 
in  the  caption.  Note  in  particular  that  a  tighter  convergence 
test  would  have  located  the  computed  minimum  closer  to  the  exact 
one.  No  studies  of  the  effects  of  variations  in  these  parameters 
were  performed  at  this  time,  since  it  was  decided  to  concentrate 
the  remaining  efforts  on  more  complicated  problems  with  program 
SLACK. 

6 . 2  Biconvex  Sandwich  Wing 

The  structure  chosen  for  the  remaining  examples  was  the 
biconvex  sandwich  wing  of  reference  9,  example  2.  A  planform 
view  of  this  wing  is  given  in  figure  7,  which  shows  the  node-point 
and  cover-panel  numbering.  The  coordinates  of  the  nodes  are  given 
in  reference  9.  To  model  the  core,  thick  shear  webs  were  placed 
along  each  chordwise  and  spanwise  line  between  node  points,  and 
axial  elements  were  placed  between  upper-surface  and  lower-surface 
node  points,  except  at  the  wing  root.  These  elements  were  used 
to  contribute  stiffness  only.  Triangular  cover  panels,  numbered 
as  indicated  in  figure  7  for  the  upper  surface,  were  used  to  model 
the  skin.  The  complete  model  had  76  nodes,  14  of  which  were 
restrained  at  the  root;  102  cover  elements;  49  shear  webs;  and  31 
axial  elements.  There  were  186  unrestrained  degrees  of  freedom, 
which  were  collapsed  to  31  transverse  degrees  of  freedom  at  the 
upper  surface.  The  design  variables  were  taken  to  be  the  thicknesses 
of  the  cover  panels.  Each  design  variable  governed  the  thickness 


of  an  upper  cover  panel  and  its  lower-surface  counterpart,  so 
there  were  a  total  of  51  design  variables,  which  could  be 
additionally  linked  in  any  combination  to  create  smaller  problems. 
Concentrated  tuning  masses  were  also  located  at  nodes  18  and  56; 
these  were  controlled  by  a  single  additional  concentrated-mass 
design  variable.  Finally,  nonactive  concentrated  masses  were 
located  at  all  the  nodes  to  represent  internal  fuel.  The  values 
used  for  these  masses  were  obtained  by  subtracting  the  active  mass 
contributed  by  the  cover  panels  from  the  total  node-point  masses 
given  in  reference  9. 

Material  properties  for  the  cover  and  web  elements  were  the 

8 

same  as  those  given  in  reference  9:  E  =  1.1308  *  10  kPa, 

p  =  4. 4244  x  lo-6  kg/cm3,  and  \T=0.3  for  the  covers;  and 

G  =  4.3492  *  107  kPa  and  ^=0.3  for  the  webs.  The  axial  elements, 

which  had  no  counterparts  in  the  model  of  reference  9,  had 
4  2 

EA  -  2.1440  *  10  kPa-m  .  These  elements  were  needed  in  the  present 
model  to  stabilize  the  shear  webs,  which  were  different  from  the 
webs  used  in  reference  9.  It  should  be  noted  that  the  webs  and 
axial  elements  were  used  to  provide  an  extremely  stiff,  weightless 
core.  As  long  as  these  elements  are  stiff  enough,  the  transverse 
mode  shapes  and  frequencies  of  the  wing  will  not  be  very  sensitive 
to  the  values  used  for  these  stiffnesses.  The  webs  were  all 
assumed  to  be  9.8  in.  thick,  and  the  cover  panels  were  sized 
according  to  the  initial  design  of  example  2  in  reference  9  - 
0.2540  cm  for  panels  1-43  and  52-94,  and  0.05080  cm  for  panels 
44-51  and  95-102.  The  initial  weight  for  this  wing  was  67,000  kg, 
including  the  fuel,  with  7,602  kg  available  for  optimization 
This  includes  350.3  kg  for  the  tuning  masses. 

In  the  examples  that  follow,  the  design  variables  were 
linked  as  indicated  in  figure  8,  so  that  the  number  of  design 
variables  was  reduced  to  seven,  including  the  concentrated-mass 
design  variable  (no.  7). 
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6.3  Optimization  with  Frequency  Constraints 

Mode  shapes  and  frequencies  were  computed  for  the  wing 

described  above,  and  the  wing  was  then  optimized  by  program 

SLACK  with  the  requirement  that  the  first  two  natural  frequencies 

be  greater  than  or  equal  to  6.20  Hz  and  18.70  Hz,  respectively. 

The  corresponding  frequencies  of  the  initial  design  were  slightly 

in  excess  of  these  values.  The  first  12  modes  of  the  initial 

design  were  used  as  fixed  modes  to  define  generalized  coordinates. 

After  34  iterations,  an  optimum  weight  (less  fuel)  of  569.9  kg 

was  computed.  An  iteration  history  of  the  variable  weight  is 

presented  in  figure  9,  and  the  optimal  values  of  the  design 

variables  are  presented  in  figure  10.  Around  25%  of  the  weight 

available  for  optimization  was  removed.  Both  frequency  constraints 

were  active  at  the  optimum,  as  would  be  expected.  CPU  time 

required  on  an  IBM  370/78  computer  was  approximately  0.1  sec  per 

iteration.  Parameters  used  in  SLACK  were  6.  =  6  =0.5,  t.  =0.9, 

t  c  x 

d  =  2.5,  K  =  0 . 1 ,  and  ev  =  0.001.  Minimum  values  of  0.02540  cm  were 
specified  for  design  variables  1-5,  while  design  variables  6  and 
7  had  minima  of  0.01016  cm  and  35.03  kg,  respectively. 

The  iteration  history  of  figure  9  demonstrates  the  usual 
rapid  reduction  in  weight,  followed  by  relatively  small  weight 
reductions  as  the  optimum  solution  is  approached. 

6.4  Optimization  with  Flutter  Constraints 

The  same  wing  was  next  analyzed  for  flutter.  The  initial 
design,  with  12  modes  again  used,  had  a  critical  flutter  speed 
of  264  m/sec  at  an  altitude  of  2,438  m,  which  was  a  matched  flutter 
point  for  Mw  =  0.8.  Doublet-lattice  aerodynamic  generalized  forces 
were  computed  with  a  code  based  on  that  used  in  reference  13. 
Polynomial  approximations  to  the  mode  shapes  were  used,  and 
computed  aerodynamic  loads  at  discrete  values  of  reduced  frequency 
were  fitted  with  polynomials  to  provide  input  for  the  flutter- 
analysis  code  and  later  on  for  program  SLACK,  where  derivatives  of 
these  loads  with  respect  to  reduced  frequency  were  required. 
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Optimization  was  originally  attempted  with  a  single  flutter 
constraint,  g  £  0  at  262  m/sec.  However,  after  six  iterations, 
another  root  became  critical  at  a  much  lower  speed,  and  it  was 
necessary  to  constrain  this  as  well.  The  branches  of  interest 
are  plotted  in  V-g  space  in  figure  11  for  the  design  obtained 
after  the  fifth  iteration.  The  mode  2  crossover  at  265  m/sec  is 
the  one  originally  constrained.  To  prevent  the  mode  4  crossover 
from  moving  to  a  lower  speed  than  this,  a  second  constraint  was 
imposed,  namely  g<_0  at  276  m/sec,  as  shown  in  figure  11.  This 
configuration  was  then  optimized,  and  the  results  are  given  in 
figures  12  and  13.  The  parameter  values  used  in  SLACK  were  the 
same  as  those  used  for  the  optimization  with  multiple  frequency 
constraints,  and  the  same  minimum  values  for  the  design  variables 
were  also  imposed. 

In  figure  12,  the  initial  weight  is  that  at  the  end  of  the 
fifth  iteration  from  the  initial  design  with  a  single  flutter 
constraint,  which  is  6,338  kg.  The  optimum,  found  after  another 
40  iterations,  is  5,647  kg,  which  is  74%  of  the  true  initial 
weight,  less  fuel,  of  7,602  kg.  The  optimal  skin  thickness 
distribution  shown  in  figure  13  is  quite  different  from  that  found 
for  multiple  frequency  constraints.  Approximately  1.2  sec  per 
iteration  were  required  by  program  SLACK  for  this  problem. 

6.5  Optimization  with  Flutter  and  Frequency  Constraints 

A  single  flutter  constraint  (g£0  at  262  m/sec)  and  a  single 
frequency  constraint  (co^^6.20  Hz)  were  then  imposed.  After  five 
iterations,  a  mode  4  crossover  at  a  speed  less  than  262  m/sec  was 
again  found.  A  flutter  analysis  of  the  fourth-iteration  design 
produced  mode  2  and  mode  4  crossovers  similar  to  those  in  figure  11, 
except  that  the  mode  4  crossover  was  at  305  m/sec.  Another  flutter 
constraint  was  added  -  g£0  at  305  m/sec  -  and  an  optimal  solution 
was  found  after  17  more  iterations.  The  results  for  this  problem 
are  given  in  figures  14  and  15.  The  initial  weight  in  figure  14 


was  that  of  the  fourth-iteration  design  -  6,466  kg.  The  final 
weight  was  virtually  identical  to  that  found  with  the  flutter 
constraints  alone  -  5,653  kg.  This  can  also  be  seen  in  figure  15, 
where  the  optimal  skin  thickness  distribution  is  very  similar  to 
that  given  in  figure  13.  The  previous  parameter  values  were  used 
here  in  SLACK,  and  approximately  1.5  sec  per  iteration  were 
required. 


7 .  CONCLUDING  REMARKS 

As  noted  in  the  discussion  of  the  Gauss-Seidel  iterative 
procedure  for  determining  the  Lagrange  multipliers,  the  strategems 
employed  in  program  GAUSS  are  most  effective  when  there  are  a 
large  number  of  constraints,  so  that  determining  the  set  of 
constraints  that  is  active  at  the  optimum  is  a  genuine  problem. 

In  particular,  this  should  be  the  case  when  stress  constraints 
are  added  to  the  behavioral  constraints  discussed  in  this  report. 
Adding  the  capability  to  compute  stresses  and  their  derivatives 
is  therefore  a  prime  requirement  in  order  to  test  program  GAUSS 
adequately. 

The  results  from  the  optimization  problems  treated  with 
program  SLACK  suggest  that  this  program  is  reliable  and  easy  to 
use.  Convergence  properties,  however,  could  be  enhanced,  probably 
by  a  better  selection  of  the  various  parameters  that  need  to  be 
supplied  by  the  user.  More  studies  with  variations  of  these 
parameters  would  be  useful.  Also,  it  is  time  to  treat  larger 
problems,  and  particularly  to  consider  problems  with  stress 
constraints,  so  that  comparisons  with  program  GAUSS  can  be 
undertaken. 
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Table  1.-  Comparison  of  exact  and  finite-difference  derivatives  for  the  biconvex 
sandwich  wing  with  two  flutter  constraints  (e. ,c0)  and  a  frequency 


Table  2.-  Variation  of  iterations  required 
for  convergence  with  various 
parameters  in  SLACK  for  test 
case  of  figure  5. 
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Figure  4.-  SLACK  iteration  history  in  design  space  for  2-DV  test 

case,  with  <5t  =  0.5,  6C  *  0.5,  t #  =  0.9, 
d  =  2.5,  K  “  2.0,  and  ev  =  0.001. 
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Figure  9.-  SLACK  weight  iteration  history  for  biconvex  wing 

with  two  frequency  constraints. 


for  biconvex  wing  with  two  frequency  constraints. 
Skin  thicknesses  in  cm;  tuning-mass  value  in  kg. 
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Figure  13.-  Optimal  values  of  design  variables  (one  surface) 
for  biconvex  wing  with  two  flutter  constraints. 
Skin  thicknesses  in  cm;  tuning-mass  value  in  kg. 


Weight  ratio,  (M^-^ 


Figure  15.-  Optimal  values  of  design  variables  (one  surface) 
for  biconvex  wing  with  two  flutter  constraints 
and  one  frequency  constraint.  Skin  thicknesses 
in  cm;  tuning-mass  value  in  kg. 
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NOMENCLATURE 


b 

c,c, 

c  ! 


C. 

1 


D  . 

3 

E 

EA 

[GA] 

[GM] , [GK] 
[GMil,  [GKJ 

k 

K 


IK.] 


M 

o 

“r 

M 


Wright  coefficient  -  see  equation  Cl) 

Preliminary  set  of  active  or  violated  constraints 

Set  of  active  or  violated  constraints  as  identified 
by  Gauss-Seidel  iteration 

Set  of  active  design  variables 

Reference  length 

Inequality  constraint  function?  jth  inequality 
constraint  function 

2 

jth  equality  constraint  function,  c^  +  z  ^ 

Redesign  factor  on  design  variables 
Limiting  value  on  |  ot— 1 1  -  see  equation  (33) 

Redesign  factor  on  slack  variables 
Young's  modulus 

Product  of  Young’s  modulus  and  cross-sectional  area 
Generalized  aerodynamics  matrix 

Generalized  mass  and  stiffness  matrices,  respectively 

Derivatives  of  generalized  mass  and  stiffness  matrices 
with  respect  to  design  variable  t^ 

Reduced  frequency,  (i>b/Vm 

Step-size  parameter  for  z^ 

Derivative  of  discrete  stiffness  matrix  with  respect 
to  design  variable  t^ 

Mass  not  available  for  optimization 

Total  mass 

Number  of  design  variables 
Mach  number 


NOMENCLATURE  (Continued) 

Number  of  constraints 

Number  of  active  design  variables 

Set  of  design  variables  predicted  to  reach  a  lower 
bound  during  an  iteration 

Set  of  all  lower-bound  passive  variables 

Set  of  design  variables  predicted  to  reach  an  upper 
bound  during  an  iteration 

Set  of  all  upper-bound  passive  variables 

ith  design  variable 

Displacment  being  constrained;  rth  displacement  in 
column  of  node-point  displacements 

Parameters  for  determining  P  and  P.  respectively 

u  * 

Convergence  function  -  see  equation  (4) 

Normalized  convergence  function,  i^/a.^ 

Free-stream  speed 
jth  slack  variable 
Step-size  parameter 

Multiplier  to  determine  a  at  each  iteration  - 
see  equation  (48) 

Functions  used  in  determining  Lagrange  multipliers  - 
see  equations  (21)  and  (22) 

Augmented  function  -  see  equation  (25) 

Step-size  parameter  used  in  determining  ]  ot— 1 1  -  see 
equation  (28) 

Constraint-change  parameter  used  in  determining 
|  ot— 1 1  -  see  equation  (30) 

Constraint  tolerance  parameter  used  to  define  A  -  see 
equation  (41)  ° 
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6 


x 


V 

V 

p 


01,(0, 


ft 

(  ) 
t  ) 


V 


des 


NONMENCLATURE  (Concluded) 

Multiplier  to  determine  fi  at  each  iteration  -  see 
equation  (49)  c 

Convergence  parameter  for  v^ 

Convergence  parameter  for  c^ 

jth  Lagrange  multiplier 

Iteration  number 

Poisson’s  ratio 

Mass  density 

Frequency;  rth  natural  frequency 
Eigenvalue  from  aeroelastic  equations 
Quantity  at  vth  iteration 
Desired  quantity 
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